24 - Time Series Regression

Learning Objectives for Today’s Class

1. Conduct tests for significance of individual regression coefficients.
   
2. Use a simple linear regression model for trend adjustment.
3. Interpret regression diagnostic plots.
   
4. Create prediction intervals for individual values of the response variable.
   
5. Discuss Exam 03 (Not in R).

Required Packages

In the code chunk below, we use the pacman package to load the required packages and install them if necessary.

if(require(pacman)==FALSE) install.packages("pacman")

## Loading required package: pacman

pacman::p_load(tidyverse, 
               magrittr, # for the two-way pipe
               fpp2, # accuracy() function, Arima(), auto.arima(),
               astsa, # package is for the required data
               plotly, # for interactive graphs
               stargazer) 

1 Recap of What we Covered Last Class

1.1 Loading the J&J Data

jj = jj # assigning the jj data (from astsa) to an object of the same name
class(jj) # to check and see if the class is a ts object

## [1] "ts"

frequency(jj) # to check the frequency --> based on the print out (=4)

## [1] 4

jj

##           Qtr1      Qtr2      Qtr3      Qtr4
## 1960  0.710000  0.630000  0.850000  0.440000
## 1961  0.610000  0.690000  0.920000  0.550000
## 1962  0.720000  0.770000  0.920000  0.600000
## 1963  0.830000  0.800000  1.000000  0.770000
## 1964  0.920000  1.000000  1.240000  1.000000
## 1965  1.160000  1.300000  1.450000  1.250000
## 1966  1.260000  1.380000  1.860000  1.560000
## 1967  1.530000  1.590000  1.830000  1.860000
## 1968  1.530000  2.070000  2.340000  2.250000
## 1969  2.160000  2.430000  2.700000  2.250000
## 1970  2.790000  3.420000  3.690000  3.600000
## 1971  3.600000  4.320000  4.320000  4.050000
## 1972  4.860000  5.040000  5.040000  4.410000
## 1973  5.580000  5.850000  6.570000  5.310000
## 1974  6.030000  6.390000  6.930000  5.850000
## 1975  6.930000  7.740000  7.830000  6.120000
## 1976  7.740000  8.910000  8.280000  6.840000
## 1977  9.540000 10.260000  9.540000  8.729999
## 1978 11.880000 12.060000 12.150000  8.910000
## 1979 14.040000 12.960000 14.850000  9.990000
## 1980 16.200000 14.670000 16.020000 11.610000

p = autoplot(jj) + theme_bw()
ggplotly(p)

Based on the plot, we can make three observations:

1. The data is not linear, which means that fitting a linear regression directly to this data is not prudent.
   
2. The variance was not constant as it increased over time (with larger values of the EPS).
   
3. The data is exhibiting a seasonal pattern (fourth quarter is consistently below the values for q3 and q1).

1.2 Transforming the J&J TS

log_jj = log(jj)
p2 = autoplot(log_jj)
ggplotly(p2) # a linear regression line is probably okay (the variance is more stable with the transformation)

1.3 Time as the Independent Variable

t = time(log_jj) # time makes a decimal date from the ts (if freq  >  1)
t

##         Qtr1    Qtr2    Qtr3    Qtr4
## 1960 1960.00 1960.25 1960.50 1960.75
## 1961 1961.00 1961.25 1961.50 1961.75
## 1962 1962.00 1962.25 1962.50 1962.75
## 1963 1963.00 1963.25 1963.50 1963.75
## 1964 1964.00 1964.25 1964.50 1964.75
## 1965 1965.00 1965.25 1965.50 1965.75
## 1966 1966.00 1966.25 1966.50 1966.75
## 1967 1967.00 1967.25 1967.50 1967.75
## 1968 1968.00 1968.25 1968.50 1968.75
## 1969 1969.00 1969.25 1969.50 1969.75
## 1970 1970.00 1970.25 1970.50 1970.75
## 1971 1971.00 1971.25 1971.50 1971.75
## 1972 1972.00 1972.25 1972.50 1972.75
## 1973 1973.00 1973.25 1973.50 1973.75
## 1974 1974.00 1974.25 1974.50 1974.75
## 1975 1975.00 1975.25 1975.50 1975.75
## 1976 1976.00 1976.25 1976.50 1976.75
## 1977 1977.00 1977.25 1977.50 1977.75
## 1978 1978.00 1978.25 1978.50 1978.75
## 1979 1979.00 1979.25 1979.50 1979.75
## 1980 1980.00 1980.25 1980.50 1980.75

1.4 Fit the Model

model = lm(log_jj ~ t) # t is the ind variable and log_jj is the response
names(model)

##  [1] "coefficients"  "residuals"     "effects"       "rank"         
##  [5] "fitted.values" "assign"        "qr"            "df.residual"  
##  [9] "xlevels"       "call"          "terms"         "model"

summary(model)

## 
## Call:
## lm(formula = log_jj ~ t)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.38309 -0.08569  0.00297  0.09984  0.38016 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -3.275e+02  5.623e+00  -58.25   <2e-16 ***
## t            1.668e-01  2.854e-03   58.45   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.1585 on 82 degrees of freedom
## Multiple R-squared:  0.9766, Adjusted R-squared:  0.9763 
## F-statistic:  3416 on 1 and 82 DF,  p-value: < 2.2e-16

Let us recap some of the outputs.

1.4.1 Regression Equation

intercept = summary(model) %>% .[['coefficients']] %>% .[1, 'Estimate']
beta1 = summary(model) %>% .[['coefficients']] %>% .[2, 'Estimate']
sigma = summary(model) %>% .[['sigma']]

\(y_t = -327.5476 + 0.1668t + \varepsilon, \, \text{where } \varepsilon\sim\mathcal N(0, ~ 0.1585^2)\)

1.4.2 Predicted/Fitted Values and Residuals

fit = model$fitted.values
# In class on Monday, I believe we examined Q3 of 1972
fit[51]

##       51 
## 1.458882

res = model$residuals 

1.4.3 Multiple R2 [NEW]

modelSummary = summary(model)
r2 = modelSummary$r.squared

\(r^2\) = 0.9765612 measures the variability of the observations (log of the EPS for J&J) around the fitted regression line.

1.4.4 Tests for Significance

We are testing the following:

\(H_0: \, \beta_1 = 0\), and our \(H_1: \, \beta_1 \ne 0\)

Given that our df = 82 and the estimated value for our \(\beta_1\), we have a t-value of 58.4506005 indicating that we are 58.4506005 standard errors above 0. Hence, the \(p\)_value is significant.

Given that \(\beta_1 \ne 0\), we can conclude that year is a significant linear predictor of the logged EPS of J&J.

1.4.5 Predicting 1981

One option is to just use the equation and plug the values for the year manually:

Q1: log of EPS for 1981 = 2.8766545; EPS = 17.7547757 Q2: log of EPS for 1981 = 2.9183537; EPS = 18.5107882

Alternatively, you can use the forecast function from the fpp2 package.

forecast(model, h=4, newdata = c(1981.0, 1981.25, 1981.5, 1981.75), level = 95)

## Warning in forecast.lm(model, h = 4, newdata = c(1981, 1981.25, 1981.5, :
## newdata column names not specified, defaulting to first variable required.

##   Point Forecast    Lo 95    Hi 95
## 1       2.876655 2.553716 3.199593
## 2       2.918354 2.595147 3.241561
## 3       2.960053 2.636572 3.283534
## 4       3.001752 2.677991 3.325513

df = data.frame(actual = log_jj, fit = model$fitted.values)
df %>% ggplot(aes(x = time(actual))) + 
  geom_line(aes(y = fit)) +
  geom_line(aes(y = actual))

## Don't know how to automatically pick scale for object of type ts. Defaulting to continuous.

2 Use a simple linear regression model for trend adjustment.

We will be using a tslm() from the fpp2 package. This allows us to automatically capture time; however, time will be coded as an integer.

modelTS = tslm(log_jj ~ trend)
summary(modelTS)

## 
## Call:
## tslm(formula = log_jj ~ trend)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.38309 -0.08569  0.00297  0.09984  0.38016 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -0.6677756  0.0349073  -19.13   <2e-16 ***
## trend        0.0416992  0.0007134   58.45   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.1585 on 82 degrees of freedom
## Multiple R-squared:  0.9766, Adjusted R-squared:  0.9763 
## F-statistic:  3416 on 1 and 82 DF,  p-value: < 2.2e-16

tsForecast = forecast(modelTS, h =8, level = 95)
autoplot(tsForecast) + autolayer(fitted(tsForecast)) + theme(legend.position = "none")

3 Interpret regression diagnostic plots.

source('resPlotTS.R') # would work assuming that you have placed the R file in the same directory

resplot(res = res, fit = fit, freq = 4)

## Loading required package: ggpubr

## 
## Attaching package: 'ggpubr'

## The following object is masked from 'package:forecast':
## 
##     gghistogram

4 Create prediction intervals for individual values of the response variable.

See the output from “Use a simple linear regression model for trend adjustment”.

5 Discuss Exam 03 (Not in R)

We have discussed until how the Limits of the ACF Plot are computed.